# 

# 

# **Metropolis-Hastings算法详解与Python代码示例**

# 

# **一、算法详解**

# 

# Metropolis-Hastings算法（简称M-H算法）是统计学和统计物理中一种重要的马尔科夫蒙特卡洛（MCMC）方法。该方法主要用于在难以直接采样的概率分布中抽取随机样本序列。其核心思想是通过构建一个马尔可夫链，使其平稳分布与目标分布相匹配，从而生成服从目标分布的样本。

# 

# 算法的名称来源于两位科学家：美国物理学家尼古拉斯·梅特罗波利斯（Nicholas Metropolis）和加拿大统计学家W·K·黑斯廷斯（W. K. Hastings）。M-H算法的基本步骤包括：

# 

# 1. **初始化**：选择一个初始状态作为马尔可夫链的起点。

# 2. **迭代过程**：

#    - **生成候选状态**：根据一个容易抽样的建议分布（proposal distribution）生成一个候选状态。

#    - **计算接受概率**：根据目标分布和建议分布计算接受候选状态的概率。

#    - **接受或拒绝**：根据接受概率决定是否接受候选状态作为下一步的状态。如果接受，则更新当前状态；如果拒绝，则保持当前状态不变。

# 3. **重复迭代**：不断重复上述过程，生成一个马尔可夫链的样本序列。

# 

# M-H算法的关键在于设计合适的建议分布和接受概率，以确保马尔可夫链的平稳分布与目标分布相匹配。此外，算法的性能还受到收敛速度和样本自相关性的影响。

# 

# **二、Python代码示例**

# 

# 下面是一个使用Python实现Metropolis-Hastings算法对一维高斯分布进行抽样的示例代码：

# 

# 

import numpy as np

import matplotlib.pyplot as plt



# 目标分布：一维高斯分布

def target_distribution(x, mean=0, std=1):

    return np.exp(-(x - mean) ** 2 / (2 * std ** 2)) / (std * np.sqrt(2 * np.pi))



# 建议分布：对称的高斯分布

def proposal_distribution(x_current, sigma=1):

    return np.random.normal(x_current, sigma)



# 接受概率

def acceptance_probability(x_current, x_candidate, target_pdf, proposal_pdf, sigma):

    A = min(1, target_pdf(x_candidate) * proposal_pdf(x_current, sigma) /

            (target_pdf(x_current) * proposal_pdf(x_candidate, sigma)))

    return A



# Metropolis-Hastings算法

def metropolis_hastings(num_samples, initial_state, target_pdf, proposal_pdf, sigma):

    samples = np.zeros(num_samples)

    current_state = initial_state

    for i in range(num_samples):

        candidate_state = proposal_distribution(current_state, sigma)

        A = acceptance_probability(current_state, candidate_state, target_pdf, proposal_pdf, sigma)

        if np.random.rand() < A:

            current_state = candidate_state

        samples[i] = current_state

    return samples



# 示例：对一维高斯分布进行抽样

num_samples = 10000

initial_state = 0

samples = metropolis_hastings(num_samples, initial_state, target_distribution, proposal_distribution, sigma=1)



# 绘制样本分布直方图

plt.hist(samples, bins=50, density=True, alpha=0.6, color='g')

# 绘制目标分布曲线

x = np.linspace(-5, 5, 100)

y = target_distribution(x)

plt.plot(x, y, 'r', linewidth=2)

plt.show()
